GLMM example 7

Preparations

Load the necessary libraries

library(car) # for regression diagnostics
library(broom) # for tidy output
library(broom.mixed) # for tidy output
library(ggfortify) # for model diagnostics
library(sjPlot) # for outputs
library(knitr) # for kable
library(effects) # for partial effects plots
library(ggeffects) # for effects plots in ggplot
library(emmeans) # for estimating marginal means
library(MASS) # for glm.nb
library(MuMIn) # for AICc
library(tidyverse) # for data wrangling
library(DHARMa) # for assessing dispersion etc
library(lme4) # for lmer
library(lmerTest) # for degrees of freedom in lmer
library(glmmTMB) # for glmmTMB
library(performance) # for diagnostic plots
library(see) # for diagnostic plots

Scenario

In an honours thesis from (1992), Mullens was investigating the ways that cane toads ( Bufo marinus ) respond to conditions of hypoxia. Toads show two different kinds of breathing patterns, lung or buccal, requiring them to be treated separately in the experiment. Her aim was to expose toads to a range of O₂ concentrations, and record their breathing patterns, including parameters such as the expired volume for individual breaths. It was desirable to have around 8 replicates to compare the responses of the two breathing types, and the complication is that animals are expensive, and different individuals are likely to have different O₂ profiles (leading to possibly reduced power). There are two main design options for this experiment;

• One animal per O₂ treatment, 8 concentrations, 2 breathing types. With 8 replicates the experiment would require 128 animals, but that this could be analysed as a completely randomised design
• One O₂ profile per animal, so that each animal would be used 8 times and only 16 animals are required (8 lung and 8 buccal breathers)

Mullens decided to use the second option so as to reduce the number of animals required (on financial and ethical grounds). By selecting this option, she did not have a set of independent measurements for each oxygen concentration, by repeated measurements on each animal across the 8 oxygen concentrations.

Toad

Format of mullens.csv data file

BREATH	TOAD	O2LEVEL	FREQBUC	SFREQBUC
lung	a	0	10.6	3.256
lung	a	5	18.8	4.336
lung	a	10	17.4	4.171
lung	a	15	16.6	4.074
...	...	...	...	...

BREATH	Categorical listing of the breathing type treatment (buccal = buccal breathing toads, lung = lung breathing toads). This is the between subjects (plots) effect and applies to the whole toads (since a single toad can only be one breathing type - either lung or buccal). Equivalent to Factor A (between plots effect) in a split-plot design
TOAD	These are the subjects (equivalent to the plots in a split-plot design: Factor B). The letters in this variable represent the labels given to each individual toad.
O2LEVEL	0 through to 50 represent the the different oxygen concentrations (0% to 50%). The different oxygen concentrations are equivalent to the within plot effects in a split-plot (Factor C).
FREQBUC	The frequency of buccal breathing - the response variable
SFREQBUC	Square root transformed frequency of buccal breathing - the response variable

Read in the data

mullens <- read_csv("../public/data/mullens.csv", trim_ws = TRUE)
mullens %>% glimpse()

## Rows: 168
## Columns: 5
## $ BREATH   <chr> "lung", "lung", "lung", "lung", "lung", "lung", "lung", "lung…
## $ TOAD     <chr> "a", "a", "a", "a", "a", "a", "a", "a", "b", "b", "b", "b", "…
## $ O2LEVEL  <dbl> 0, 5, 10, 15, 20, 30, 40, 50, 0, 5, 10, 15, 20, 30, 40, 50, 0…
## $ FREQBUC  <dbl> 10.6, 18.8, 17.4, 16.6, 9.4, 11.4, 2.8, 4.4, 21.6, 17.4, 22.4…
## $ SFREQBUC <dbl> 3.255764, 4.335897, 4.171331, 4.074310, 3.065942, 3.376389, 1…

Exploratory data analysis

Model formula: \[ y_i \sim{} \mathcal{Pois}(\lambda_i)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]

where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of copper, distance and their interaction on the number of number of worms. Area of the place segment was also incorporated as an offset. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual plates.

In this example, the individual TOADs are the random blocks. Each toad can only be of one BREATH type (they are either predominantly buccal breathers or predominantly lung breathers). Hence, BREATH is a between TOAD (block) effect. The frequency of buccal breathing of each TOAD was measured under eight different oxygen levels and thus, these represent the within block effect, as will the interaction between breathing type and oxygen level.

The response in this case is a little tricky since it is a proportion, but without the full original measurements. For example a frequency of 50% would indicate that half of the breaths taken by the toad during the monitoring phase were buccal breaths. Ideally, it would be good to have the actual counts (both the number of buccal breaths and the total number of breaths). That way, we could model the data against a binomial distribution.

As it is, we only have the proportion (as a percentage). Although we could model this against a beta distribution, this could be complicated if there are proportions of either 0 or 1 (100).

To help guide us through this and the other typical model assumptions, lets start by graphing the frequency of buccal breathing against breathing type and oxygen level. However, before we do, we need to make sure that all categorical variables are declared as factors (including the random effect). If we intend to model against a beta distribution, we will need a version of the response that is represented on a scale between 0 and 1 (but not include 0 or 1).

mullens <- mullens %>%
  mutate(
    BREATH = factor(BREATH),
    TOAD = factor(TOAD),
    pBUC = FREQBUC / 100,
    pzBUC = ifelse(pBUC == 0, 0.01, pBUC)
  )

So starting with the raw response.

ggplot(mullens, aes(y = FREQBUC, x = factor(O2LEVEL), color = BREATH)) +
  geom_boxplot()

Conclusions:

• there is a very clear relationship between mean and variance (boxplots that are higher up the y-axis are taller).
• it might be possible to address this via a logarithmic or root transformation, however this is the least favourable option (particularly root transformations due to issues with back-transformations).

If we intend to use a beta distribution, we could repeat the above with the scaled response. Since the rescaling maintains the same ranking and relative spacing, the plot will look the same as above, just with a different y-axis. However, our interpretation will change. Under a beta distribution (with logit link), we expect that the variance and mean will be related. We expect that the distributions will become more varied and more symmetrical as the expected mean shifts towards 0.5. Distributions approaching 0 and 1 will be more asymmetrical and smaller. This indeed does seem to be the case here.

Now lets explore the oxygen trends (separately for each breathing type).

ggplot(mullens, aes(y = pzBUC, x = O2LEVEL, color = BREATH)) +
  geom_smooth() +
  geom_point()

Conclusions:

• although linearity seems reasonable for the buccal breathers, there is definitely evidence of non linearity in the lung breathers.
• it might be interesting to fit polynomial trends for oxygen concentration.

ggplot(mullens, aes(y = pzBUC, x = O2LEVEL, color = BREATH)) +
  geom_smooth() +
  geom_point() +
  facet_wrap(~ BREATH + TOAD, scales = "free")

## facet_grid(TOAD~BREATH)

Conclusions:

• it does appear that different toads (within the same breathing type) have very different levels of buccal breathing
• the individual toads also have varying responses to changes in oxygen concentration
• it might be useful to explore random intercept/slope models.

Fit the model

sqrt (Gaussian) with lmer

random structure

mullens.lmer1a <- lmer(SFREQBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD),
  data = mullens,
  REML = TRUE
)
mullens.lmer1b <- lmer(SFREQBUC ~ BREATH * poly(O2LEVEL, 3) + (poly(O2LEVEL, 3) | TOAD),
  data = mullens,
  REML = TRUE
)
# OR
mullens.lmer1b <- update(mullens.lmer1a, . ~ . - (1 | TOAD) + (poly(O2LEVEL, 3) | TOAD))
mullens.allFit <- allFit(mullens.lmer1b)

## bobyqa : [OK]
## Nelder_Mead : [OK]
## nlminbwrap : [OK]
## nmkbw : [OK]
## optimx.L-BFGS-B : [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
## nloptwrap.NLOPT_LN_BOBYQA : [OK]

mullens.allFit

## original model:
## SFREQBUC ~ BREATH + poly(O2LEVEL, 3) + (poly(O2LEVEL, 3) | TOAD) + BREATH:pol... 
## data:  mullens 
## optimizers (7): bobyqa, Nelder_Mead, nlminbwrap, nmkbw, optimx.L-BFGS-B,nloptwrap.NLOPT_LN_N...
## differences in negative log-likelihoods:
## max= 0.987 ; std dev= 0.364

## Check which of the models are considered valid (OK)
is.OK <- sapply(mullens.allFit, is, "merMod")
is.OK

##                        bobyqa                   Nelder_Mead 
##                          TRUE                          TRUE 
##                    nlminbwrap                         nmkbw 
##                          TRUE                          TRUE 
##               optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD 
##                          TRUE                          TRUE 
##     nloptwrap.NLOPT_LN_BOBYQA 
##                          TRUE

AICc(mullens.lmer1a, mullens.lmer1b)

## mullens.lmer1c <- update(mullens.lmer1a,  ~. -(1|TOAD) + (poly(O2LEVEL, 3)||TOAD))
## anova(mullens.lmer1c)
## mullens.lmer1d <- update(mullens.lmer1a,  ~. + (0+poly(O2LEVEL, 3)||TOAD))
## anova(mullens.lmer1d)

Conclusions:

• the more complex random intercept/slope model failed to converge.
• essentially this is because the negative log-likelihood surface does not have a discernible trough (minimum). Instead it is either flat or monotonically declining. Hence the optimiser has not identified a minimum (other than at the end of the surface).
• this is usually because the model is over-parameterised (too complex for the amount of data available).
• it can also be due to random effects (or dispersion or zero-inflation) that are estimated to be zero.
• for Bayesian models we could tighten the priors and try refitting, but for Frequentest models, there is no alternative but to proceed with the more simple model.

fixed structure

mullens.lmer1c <- lmer(SFREQBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD),
  data = mullens,
  REML = FALSE
)
mullens.lmer1d <- update(mullens.lmer1c, . ~ . - BREATH:poly(O2LEVEL, 3))
AICc(mullens.lmer1c, mullens.lmer1d)

Conclusions:

• proceed with multiplicative model

sqrt (Gaussian)

random structure

mullens.glmmTMB1a <- glmmTMB(SFREQBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD),
  data = mullens,
  family = gaussian(),
  REML = TRUE
)
mullens.glmmTMB1b <- glmmTMB(SFREQBUC ~ BREATH * poly(O2LEVEL, 3) + (poly(O2LEVEL, 3) | TOAD),
  data = mullens,
  family = gaussian(),
  REML = TRUE,
  control = glmmTMBControl(
    optimizer = optim,
    optArgs = list(method = "BFGS")
  )
)
## OR
mullens.glmmTMB1b <- update(mullens.glmmTMB1a, ~ . - (1 | TOAD) + (poly(O2LEVEL, 3) | TOAD),
  control = glmmTMBControl(
    optimizer = optim,
    optArgs = list(method = "BFGS")
  )
)
AICc(mullens.glmmTMB1a, mullens.glmmTMB1b)

Conclusions:

• the more complex random intercept/slope model failed to converge.
• essentially this is because the negative log-likelihood surface does not have a discernible trough (minimum). Instead it is either flat or monotonically declining. Hence the optimiser has not identified a minimum (other than at the end of the surface).
• this is usually because the model is over-parameterised (too complex for the amount of data available).
• it can also be due to random effects (or dispersion or zero-inflation) that are estimated to be zero.
• for Bayesian models we could tighten the priors and try refitting, but for Frequentest models, there is no alternative but to proceed with the more simple model.

fixed structure

## mullens.glmmTMB1c <- glmmTMB(SFREQBUC ~ BREATH*poly(O2LEVEL, 3) + (1|TOAD), data=mullens,
##                             family=gaussian(),  REML=FALSE)
## mullens.glmmTMB1d = update(mullens.glmmTMB1a,  .~.-BREATH:poly(O2LEVEL,  3))
## AICc(mullens.glmmTMB1a,  mullens.glmmTMB1b)

Conclusions:

• proceed with multiplicative model

beta

random structure

mullens.glmmTMB2a <- glmmTMB(pzBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD),
  data = mullens,
  family = beta_family(link = "logit"),
  REML = TRUE
)
mullens.glmmTMB2b <- glmmTMB(pzBUC ~ BREATH * poly(O2LEVEL, 3) + (poly(O2LEVEL, 3) | TOAD),
  data = mullens,
  family = beta_family(link = "logit"),
  REML = TRUE,
  control = glmmTMBControl(
    optimizer = optim,
    optArgs = list(method = "BFGS")
  )
)
## OR
mullens.glmmTMB2b <- update(mullens.glmmTMB2a, ~ . - (1 | TOAD) + (poly(O2LEVEL, 3) | TOAD),
  control = glmmTMBControl(
    optimizer = optim,
    optArgs = list(method = "BFGS")
  )
)
AICc(mullens.glmmTMB2a, mullens.glmmTMB2b)

Conclusions:

• proceed with the more complex random intercept/random slope model

fixed structure

## mullens.glmmTMB2a = glmmTMB(pzBUC ~ BREATH*poly(O2LEVEL, 3) + (1|TOAD), data=mullens,
##                             family=beta_family(link = "logit"),  REML=FALSE)
## mullens.glmmTMB2b = update(mullens.glmmTMB2a,  .~.-BREATH:poly(O2LEVEL,  3))
## AICc(mullens.glmmTMB2a,  mullens.glmmTMB2b)

Conclusions:

• proceed with multiplicative model (it has a more negative AICc)

Model validation

sqrt (Gaussian) with lmer

plot_model

plot_model(mullens.lmer1a, type = "diag")[-2] %>% plot_grid()

Performance model checking

mullens.lmer1a %>% performance::check_model()

DHARMa residuals

mullens.resid <- mullens.lmer1a %>% simulateResiduals(plot = TRUE)

sqrt (Gaussian)

plot_model

plot_model(mullens.glmmTMB1b, type = "diag")[-2] %>% plot_grid()

Performance model checking

mullens.glmmTMB1b %>% performance::check_model()

DHARMa residuals

mullens.resid <- mullens.glmmTMB1b %>% simulateResiduals(plot = TRUE)

beta

plot_model

plot_model(mullens.glmmTMB2b, type = "diag")[-2] %>% plot_grid()

Performance model checking

mullens.glmmTMB2b %>% performance::check_model()

DHARMa residuals

## mullens.resid <- mullens.glmmTMB2a %>% simulateResiduals(plot=TRUE)
mullens.resid <- mullens.glmmTMB2b %>% simulateResiduals(plot = TRUE)

mullens.glmmTMB2b %>%
  ggemmeans(~ O2LEVEL | BREATH) %>%
  plot(add.data = TRUE)

Model investigation / hypothesis testing

sqrt (Gaussian) with lmer

summary

mullens.lmer1a %>% summary()

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SFREQBUC ~ BREATH * poly(O2LEVEL, 3) + (1 | TOAD)
##    Data: mullens
## 
## REML criterion at convergence: 453.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3421 -0.5127  0.0282  0.4435  3.4313 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TOAD     (Intercept) 0.7733   0.8794  
##  Residual             0.7416   0.8612  
## Number of obs: 168, groups:  TOAD, 21
## 
## Fixed effects:
##                               Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)                    3.77246    0.25810  19.00000  14.616 8.68e-12
## BREATHlung                    -1.00381    0.41817  19.00000  -2.400   0.0268
## poly(O2LEVEL, 3)1            -10.76355    1.09451 141.00000  -9.834  < 2e-16
## poly(O2LEVEL, 3)2              1.31085    1.09451 141.00000   1.198   0.2331
## poly(O2LEVEL, 3)3             -0.01429    1.09451 141.00000  -0.013   0.9896
## BREATHlung:poly(O2LEVEL, 3)1  13.03414    1.77330 141.00000   7.350 1.47e-11
## BREATHlung:poly(O2LEVEL, 3)2  -7.22946    1.77330 141.00000  -4.077 7.59e-05
## BREATHlung:poly(O2LEVEL, 3)3   2.75036    1.77330 141.00000   1.551   0.1231
##                                 
## (Intercept)                  ***
## BREATHlung                   *  
## poly(O2LEVEL, 3)1            ***
## poly(O2LEVEL, 3)2               
## poly(O2LEVEL, 3)3               
## BREATHlung:poly(O2LEVEL, 3)1 ***
## BREATHlung:poly(O2LEVEL, 3)2 ***
## BREATHlung:poly(O2LEVEL, 3)3    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                     (Intr) BREATH p(O2LEVEL,3)1 p(O2LEVEL,3)2 p(O2LEVEL,3)3
## BREATHlung          -0.617                                                 
## p(O2LEVEL,3)1        0.000  0.000                                          
## p(O2LEVEL,3)2        0.000  0.000  0.000                                   
## p(O2LEVEL,3)3        0.000  0.000  0.000         0.000                     
## BREATH:(O2LEVEL,3)1  0.000  0.000 -0.617         0.000         0.000       
## BREATH:(O2LEVEL,3)2  0.000  0.000  0.000        -0.617         0.000       
## BREATH:(O2LEVEL,3)3  0.000  0.000  0.000         0.000        -0.617       
##                     BREATH:(O2LEVEL,3)1 BREATH:(O2LEVEL,3)2
## BREATHlung                                                 
## p(O2LEVEL,3)1                                              
## p(O2LEVEL,3)2                                              
## p(O2LEVEL,3)3                                              
## BREATH:(O2LEVEL,3)1                                        
## BREATH:(O2LEVEL,3)2  0.000                                 
## BREATH:(O2LEVEL,3)3  0.000               0.000

• The parameters are presented on a root scale from which it can be dangerous to back-transform. Hence the parameters themselves have little value (other which were not found to be significant).
• this might imply that the lung breathing toads have a quadratic trend, but this is not formally tested here.
• there is approximately the same amount of variation between toads (within breathing types) as there is within toads (for a given oxygen concentration).

tidy

mullens.lmer1a %>% tidy(conf.int = TRUE)

mullens.lmer1a %>%
  tidy(conf.int = TRUE) %>%
  kable()

effect	group	term	estimate	std.error	statistic	df	p.value	conf.low	conf.high
fixed	NA	(Intercept)	3.7724598	0.258102	14.6161608	19	0.0000000	3.2322461	4.3126734
fixed	NA	BREATHlung	-1.0038085	0.418173	-2.4004623	19	0.0267793	-1.8790545	-0.1285624
fixed	NA	poly(O2LEVEL, 3)1	-10.7635498	1.094506	-9.8341613	141	0.0000000	-12.9273134	-8.5997862
fixed	NA	poly(O2LEVEL, 3)2	1.3108539	1.094506	1.1976670	141	0.2330562	-0.8529097	3.4746175
fixed	NA	poly(O2LEVEL, 3)3	-0.0142947	1.094506	-0.0130604	141	0.9895981	-2.1780583	2.1494690
fixed	NA	BREATHlung:poly(O2LEVEL, 3)1	13.0341442	1.773303	7.3502087	141	0.0000000	9.5284464	16.5398419
fixed	NA	BREATHlung:poly(O2LEVEL, 3)2	-7.2294601	1.773303	-4.0768338	141	0.0000759	-10.7351578	-3.7237623
fixed	NA	BREATHlung:poly(O2LEVEL, 3)3	2.7503615	1.773303	1.5509826	141	0.1231475	-0.7553363	6.2560593
ran_pars	TOAD	sd__(Intercept)	0.8793850	NA	NA	NA	NA	NA	NA
ran_pars	Residual	sd__Observation	0.8611528	NA	NA	NA	NA	NA	NA

Conclusions:

• The parameters are presented on a root scale from which it can be dangerous to back-transform. Hence the parameters themselves have little value (other which were not found to be significant).
• this might imply that the lung breathing toads have a quadratic trend, but this is not formally tested here.
• there is approximately the same amount of variation between toads (within breathing types) as there is within toads (for a given oxygen concentration).

tab_model

# warning this is only appropriate for html output
mullens.lmer1a %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)

	SFREQBUC
Predictors	Estimates	std. Error	CI	p
(Intercept)	3.77	0.26	3.26 – 4.28	<0.001
BREATH [lung]	-1.00	0.42	-1.83 – -0.18	0.018
O2LEVEL [1st degree]	-10.76	1.09	-12.93 – -8.60	<0.001
O2LEVEL [2nd degree]	1.31	1.09	-0.85 – 3.47	0.233
O2LEVEL [3rd degree]	-0.01	1.09	-2.18 – 2.15	0.990
BREATH [lung] * O2LEVEL [1st degree]	13.03	1.77	9.53 – 16.54	<0.001
BREATH [lung] * O2LEVEL [2nd degree]	-7.23	1.77	-10.73 – -3.73	<0.001
BREATH [lung] * O2LEVEL [3rd degree]	2.75	1.77	-0.75 – 6.25	0.123
Random Effects
σ2	0.74
τ00 TOAD	0.77
ICC	0.51
N TOAD	21
Observations	168
Marginal R2 / Conditional R2	0.341 / 0.677
AIC	473.152

anova

# with Satterthwaite degrees of freedom calculations
mullens.lmer1a %>% anova()

# OR with Keyward-Roger degrees of freedom calculations
mullens.lmer1a %>% anova(ddf = "Kenward-Roger")

sqrt (Gaussian)

summary

mullens.glmmTMB1b %>% summary()

##  Family: gaussian  ( identity )
## Formula:          
## SFREQBUC ~ BREATH + poly(O2LEVEL, 3) + (poly(O2LEVEL, 3) | TOAD) +  
##     BREATH:poly(O2LEVEL, 3)
## Data: mullens
## 
##      AIC      BIC   logLik deviance df.resid 
##    452.6    511.9   -207.3    414.6      157 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name              Variance Std.Dev. Corr              
##  TOAD     (Intercept)        0.8127  0.9015                     
##           poly(O2LEVEL, 3)1 35.4538  5.9543   -0.11             
##           poly(O2LEVEL, 3)2 10.3349  3.2148   -0.30 -0.91       
##           poly(O2LEVEL, 3)3  2.3969  1.5482    0.13 -0.50  0.35 
##  Residual                    0.4320  0.6573                     
## Number of obs: 168, groups:  TOAD, 21
## 
## Dispersion estimate for gaussian family (sigma^2): 0.432 
## 
## Conditional model:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                    3.77246    0.25821  14.610  < 2e-16 ***
## BREATHlung                    -1.00381    0.41835  -2.399  0.01642 *  
## poly(O2LEVEL, 3)1            -10.76355    1.85070  -5.816 6.03e-09 ***
## poly(O2LEVEL, 3)2              1.31085    1.22182   1.073  0.28333    
## poly(O2LEVEL, 3)3             -0.01429    0.93928  -0.015  0.98786    
## BREATHlung:poly(O2LEVEL, 3)1  13.03414    2.99847   4.347 1.38e-05 ***
## BREATHlung:poly(O2LEVEL, 3)2  -7.22946    1.97958  -3.652  0.00026 ***
## BREATHlung:poly(O2LEVEL, 3)3   2.75036    1.52180   1.807  0.07071 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• The parameters are presented on a root scale from which it can be dangerous to back-transform. Hence the parameters themselves have little value (other which were not found to be significant).
• this might imply that the lung breathing toads have a quadratic trend, but this is not formally tested here.
• there is approximately the same amount of variation between toads (within breathing types) as there is within toads (for a given oxygen concentration).

tidy

mullens.glmmTMB1b %>% tidy()

mullens.glmmTMB1b %>% tidy(effects = "fixed", conf.int = TRUE)

mullens.glmmTMB1b %>%
  tidy(effects = "fixed", conf.int = TRUE) %>%
  kable()

effect	component	term	estimate	std.error	statistic	p.value	conf.low	conf.high
fixed	cond	(Intercept)	3.7724598	0.2582114	14.6099634	0.0000000	3.2663746	4.2785449
fixed	cond	BREATHlung	-1.0038085	0.4183504	-2.3994445	0.0164200	-1.8237601	-0.1838568
fixed	cond	poly(O2LEVEL, 3)1	-10.7635498	1.8506966	-5.8159452	0.0000000	-14.3908484	-7.1362512
fixed	cond	poly(O2LEVEL, 3)2	1.3108539	1.2218219	1.0728682	0.2833303	-1.0838731	3.7055808
fixed	cond	poly(O2LEVEL, 3)3	-0.0142947	0.9392759	-0.0152188	0.9878576	-1.8552416	1.8266523
fixed	cond	BREATHlung:poly(O2LEVEL, 3)1	13.0341442	2.9984711	4.3469300	0.0000138	7.1572488	18.9110396
fixed	cond	BREATHlung:poly(O2LEVEL, 3)2	-7.2294601	1.9795777	-3.6520213	0.0002602	-11.1093612	-3.3495590
fixed	cond	BREATHlung:poly(O2LEVEL, 3)3	2.7503615	1.5218009	1.8073071	0.0707144	-0.2323134	5.7330364

Conclusions:

• The parameters are presented on a root scale from which it can be dangerous to back-transform. Hence the parameters themselves have little value (other which were not found to be significant).
• this might imply that the lung breathing toads have a quadratic trend, but this is not formally tested here.
• there is approximately the same amount of variation between toads (within breathing types) as there is within toads (for a given oxygen concentration).

tab_model

# warning this is only appropriate for html output
mullens.glmmTMB1b %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)

	SFREQBUC
Predictors	Estimates	std. Error	CI	p
(Intercept)	3.77	0.26	3.27 – 4.28	<0.001
BREATH [lung]	-1.00	0.42	-1.82 – -0.18	0.016
O2LEVEL [1st degree]	-10.76	1.85	-14.39 – -7.14	<0.001
O2LEVEL [2nd degree]	1.31	1.22	-1.08 – 3.71	0.283
O2LEVEL [3rd degree]	-0.01	0.94	-1.86 – 1.83	0.988
BREATH [lung] * O2LEVEL [1st degree]	13.03	3.00	7.16 – 18.91	<0.001
BREATH [lung] * O2LEVEL [2nd degree]	-7.23	1.98	-11.11 – -3.35	<0.001
BREATH [lung] * O2LEVEL [3rd degree]	2.75	1.52	-0.23 – 5.73	0.071
Random Effects
σ2	0.43
τ00 TOAD	0.81
τ11 TOAD.poly(O2LEVEL, 3)1	35.45
τ11 TOAD.poly(O2LEVEL, 3)2	10.33
τ11 TOAD.poly(O2LEVEL, 3)3	2.40
ρ01	-0.11
	-0.30
	0.13
ICC	0.72
N TOAD	21
Observations	168
Marginal R2 / Conditional R2	0.338 / 0.813
AIC	452.576

beta

summary

mullens.glmmTMB2b %>% summary()

##  Family: beta  ( logit )
## Formula:          
## pzBUC ~ BREATH + poly(O2LEVEL, 3) + (poly(O2LEVEL, 3) | TOAD) +  
##     BREATH:poly(O2LEVEL, 3)
## Data: mullens
## 
##      AIC      BIC   logLik deviance df.resid 
##   -493.0   -433.6    265.5   -531.0      157 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name              Variance Std.Dev. Corr              
##  TOAD   (Intercept)        0.3140  0.5603                     
##         poly(O2LEVEL, 3)1 11.9893  3.4626    0.11             
##         poly(O2LEVEL, 3)2  2.3761  1.5415   -0.26 -0.99       
##         poly(O2LEVEL, 3)3  0.8209  0.9060   -0.34 -0.85  0.90 
## Number of obs: 168, groups:  TOAD, 21
## 
## Dispersion parameter for beta family (): 52.4 
## 
## Conditional model:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   -1.8176     0.1610 -11.291  < 2e-16 ***
## BREATHlung                    -0.5658     0.2622  -2.158 0.030912 *  
## poly(O2LEVEL, 3)1             -6.7617     1.1184  -6.046 1.49e-09 ***
## poly(O2LEVEL, 3)2              0.5762     0.6893   0.836 0.403184    
## poly(O2LEVEL, 3)3              0.1217     0.5932   0.205 0.837509    
## BREATHlung:poly(O2LEVEL, 3)1   7.5618     1.8474   4.093 4.26e-05 ***
## BREATHlung:poly(O2LEVEL, 3)2  -4.0459     1.1650  -3.473 0.000515 ***
## BREATHlung:poly(O2LEVEL, 3)3   1.8213     1.0269   1.774 0.076133 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Conclusions:

• these parameter estimates are on a logit (log odds/odds ratio scale)
• the estimated proportion of buccal breathing in the buccal breathing toads at an oxygen concentration of 0 is -1.82 (on the logit scale). This equates to an odds of breaths being buccal (when the toad is a buccal breather in 0% oxygen) of 0.16:1 - in other words, the odds of breaths not being buccal are 6.16:1. Hence buccal breathing toads (in 0% oxygen) are 6.16 times more likely to lung breath than buccal breath. Yet another way to express this is that buccal breathing toads in 0% oxygen buccal breath 13.97% of the time.
• lung breathing toads in 0% oxygen buccal breath -0.57 fewer times than buccal breathing toads. This equates to a 1.76 fold decline in the rate of buccal breathing between buccal and lung breathing toads at 0% oxygen.
• the estimated parameters for the polynomial terms (including the intercepts) are of little importance. However, the associated inferences are interesting
• there is evidence of a significant, negative linear trend with oxygen concentration for the buccal breathing toads.
• the degree of linearity in the lung breathing toads differs significantly from that of the buccal breathing toads.
• the degree of the quadratic component for the lung breathing toads differs from that of the buccal breathing toads (which were not found to be significant).
• this might imply that the lung breathing toads have a quadratic trend, but this is not formally tested here.
• the amount of variation between toads (within the same breathing type) is nearly twice that of the variation within toads ( 3.07 vs 0.16 )
• within toads there is a relatively high amount of variability in the linearity of the trends (compared to the variability in quadratic, cubic and residuals).

tidy

mullens.glmmTMB2b %>% tidy(conf.int = TRUE)

mullens.glmmTMB2b %>% tidy(conf.int = TRUE, exponentiate = TRUE)

mullens.glmmTMB2b %>%
  tidy(conf.int = TRUE, exponentiate = TRUE) %>%
  kable()

effect	component	group	term	estimate	std.error	statistic	p.value	conf.low	conf.high
fixed	cond	NA	(Intercept)	0.1624188	0.0261464	-11.2906065	0.0000000	0.1184701	2.226711e-01
fixed	cond	NA	BREATHlung	0.5678989	0.1488849	-2.1582031	0.0309120	0.3397138	9.493555e-01
fixed	cond	NA	poly(O2LEVEL, 3)1	0.0011573	0.0012944	-6.0456042	0.0000000	0.0001292	1.036240e-02
fixed	cond	NA	poly(O2LEVEL, 3)2	1.7793030	1.2264786	0.8359485	0.4031839	0.4608019	6.870456e+00
fixed	cond	NA	poly(O2LEVEL, 3)3	1.1293592	0.6699145	0.2050814	0.8375085	0.3531144	3.612008e+00
fixed	cond	NA	BREATHlung:poly(O2LEVEL, 3)1	1923.3855226	3553.3226762	4.0931655	0.0000426	51.4671290	7.187912e+04
fixed	cond	NA	BREATHlung:poly(O2LEVEL, 3)2	0.0174948	0.0203821	-3.4727187	0.0005152	0.0017833	1.716295e-01
fixed	cond	NA	BREATHlung:poly(O2LEVEL, 3)3	6.1799931	6.3463423	1.7735771	0.0761331	0.8257990	4.624892e+01
ran_pars	cond	TOAD	sd__(Intercept)	0.5603453	NA	NA	NA	0.3974891	7.899258e-01
ran_pars	cond	TOAD	sd__poly(O2LEVEL, 3)1	3.4625600	NA	NA	NA	2.2148406	5.413176e+00
ran_pars	cond	TOAD	sd__poly(O2LEVEL, 3)2	1.5414544	NA	NA	NA	0.7802785	3.045171e+00
ran_pars	cond	TOAD	sd__poly(O2LEVEL, 3)3	0.9060194	NA	NA	NA	0.1103516	7.438691e+00
ran_pars	cond	TOAD	cor__(Intercept).poly(O2LEVEL, 3)1	0.1053372	NA	NA	NA	-0.4110310	5.524276e-01
ran_pars	cond	TOAD	cor__(Intercept).poly(O2LEVEL, 3)2	-0.2594979	NA	NA	NA	-0.2619112	2.622242e-01
ran_pars	cond	TOAD	cor__(Intercept).poly(O2LEVEL, 3)3	-0.3406095	NA	NA	NA	-0.3386097	3.304468e-01
ran_pars	cond	TOAD	cor__poly(O2LEVEL, 3)1.poly(O2LEVEL, 3)2	-0.9862934	NA	NA	NA	-0.7721392	9.492451e-01
ran_pars	cond	TOAD	cor__poly(O2LEVEL, 3)1.poly(O2LEVEL, 3)3	-0.8485818	NA	NA	NA	-0.5798108	8.231511e-01
ran_pars	cond	TOAD	cor__poly(O2LEVEL, 3)2.poly(O2LEVEL, 3)3	0.9000072	NA	NA	NA	0.8304543	8.480642e-01

• these parameter estimates are on a logit (log odds/odds ratio scale)
• the estimated proportion of buccal breathing in the buccal breathing toads at an oxygen concentration of 0 is -1.82 (on the logit scale). This equates to an odds of breaths being buccal (when the toad is a buccal breather in 0% oxygen) of 0.16:1 - in other words, the odds of breaths not being buccal are 6.16:1. Hence buccal breathing toads (in 0% oxygen) are 6.16 times more likely to lung breath than buccal breath. Yet another way to express this is that buccal breathing toads in 0% oxygen buccal breath 13.97% of the time.
• lung breathing toads in 0% oxygen buccal breath -0.57 fewer times than buccal breathing toads. This equates to a 1.76 fold decline in the rate of buccal breathing between buccal and lung breathing toads at 0% oxygen.
• the estimated parameters for the polynomial terms (including the intercepts) are of little importance. However, the associated inferences are interesting
• there is evidence of a significant, negative linear trend with oxygen concentration for the buccal breathing toads.
• the degree of linearity in the lung breathing toads differs significantly from that of the buccal breathing toads.
• the degree of the quadratic component for the lung breathing toads differs from that of the buccal breathing toads (which were not found to be significant).
• this might imply that the lung breathing toads have a quadratic trend, but this is not formally tested here.
• the amount of variation between toads (within the same breathing type) is nearly twice that of the variation within toads ( 3.07 vs 0.16 )
• within toads there is a relatively high amount of variability in the linearity of the trends (compared to the variability in quadratic, cubic and residuals).

tab_model

# warning this is only appropriate for html output
mullens.glmmTMB2b %>% sjPlot::tab_model(show.se = TRUE, show.aic = TRUE)

	pzBUC
Predictors	Estimates	std. Error	CI	p
(Intercept)	0.16	0.03	0.12 – 0.22	<0.001
BREATH [lung]	0.57	0.15	0.34 – 0.95	0.031
O2LEVEL [1st degree]	0.00	0.00	0.00 – 0.01	<0.001
O2LEVEL [2nd degree]	1.78	1.23	0.46 – 6.87	0.403
O2LEVEL [3rd degree]	1.13	0.67	0.35 – 3.61	0.838
BREATH [lung] * O2LEVEL [1st degree]	1923.39	3553.32	51.47 – 71879.12	<0.001
BREATH [lung] * O2LEVEL [2nd degree]	0.02	0.02	0.00 – 0.17	0.001
BREATH [lung] * O2LEVEL [3rd degree]	6.18	6.35	0.83 – 46.25	0.076
Random Effects
σ2	0.16
τ00 TOAD	0.31
τ11 TOAD.poly(O2LEVEL, 3)1	11.99
τ11 TOAD.poly(O2LEVEL, 3)2	2.38
τ11 TOAD.poly(O2LEVEL, 3)3	0.82
ρ01	0.11
	-0.26
	-0.34
ICC	0.71
N TOAD	21
Observations	168
Marginal R2 / Conditional R2	0.334 / 0.809
AIC	-492.981

Further analyses

The results so far suggest that the relationship between buccal breathing frequency and oxygen concentration is different for buccal breathers compared to lung breathers. We also know that there is strong evidence for a linear trend for the buccal breathers. It would therefore be interesting to specifically test whether there is evidence for either a linear or quadratic trend for the lung breathers.

sqrt (Gaussian) with lmer

mullens.lmer1a %>%
  emtrends(specs = "BREATH", var = "O2LEVEL", max.degree = 3) %>%
  summary(infer = TRUE)

## OR
mullens.lmer1a %>% emtrends(specs = "BREATH", var = "O2LEVEL", max.degree = 3, infer = c(TRUE, TRUE))

## degree = linear:
##  BREATH O2LEVEL.trend       SE  df  lower.CL  upper.CL t.ratio p.value
##  buccal     -5.38e-02 1.33e-02 141 -8.01e-02 -2.75e-02  -4.048  0.0001
##  lung       -4.88e-03 1.69e-02 141 -3.84e-02  2.86e-02  -0.288  0.7740
## 
## degree = quadratic:
##  BREATH O2LEVEL.trend       SE  df  lower.CL  upper.CL t.ratio p.value
##  buccal      4.45e-04 4.57e-04 141 -4.59e-04  1.35e-03   0.972  0.3326
##  lung       -2.67e-03 5.83e-04 141 -3.82e-03 -1.52e-03  -4.579  <.0001
## 
## degree = cubic:
##  BREATH O2LEVEL.trend       SE  df  lower.CL  upper.CL t.ratio p.value
##  buccal     -3.50e-07 2.69e-05 141 -5.36e-05  5.29e-05  -0.013  0.9896
##  lung        6.73e-05 3.43e-05 141 -5.50e-07  1.35e-04   1.961  0.0518
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

newdata <- with(mullens, list(
  O2LEVEL = modelr::seq_range(O2LEVEL, n = 100),
  BREATH = levels(BREATH)
))
mullens.grid <- mullens.lmer1a %>%
  emmeans(~ O2LEVEL | BREATH, at = newdata) %>%
  as.data.frame()
mullens.grid %>%
  group_by(BREATH) %>%
  summarise(value = O2LEVEL[which.max(emmean)])

Conclusions:

• the rate of buccal breathing by buccal breathing toads declines linearly with increasing oxygen concentration
• the rate of buccal breathing by lung breathing toads as oxygen concentration increases follows a quadratic trend. Initially, as oxygen concentrations increases, the rate of buccal breathing increases up to an oxygen concentration of approximately 20 percent, after which it declines.

sqrt (Gaussian)

mullens.glmmTMB1b %>%
  emtrends(specs = "BREATH", var = "O2LEVEL", max.degree = 3) %>%
  summary(infer = TRUE)

## OR
emtrends(mullens.glmmTMB1b, specs = "BREATH", var = "O2LEVEL", max.degree = 3, infer = c(TRUE, TRUE))

## degree = linear:
##  BREATH O2LEVEL.trend       SE  df  lower.CL  upper.CL t.ratio p.value
##  buccal     -5.38e-02 1.63e-02 157 -8.60e-02 -2.17e-02  -3.306  0.0012
##  lung       -4.88e-03 2.07e-02 157 -4.59e-02  3.61e-02  -0.235  0.8145
## 
## degree = quadratic:
##  BREATH O2LEVEL.trend       SE  df  lower.CL  upper.CL t.ratio p.value
##  buccal      4.45e-04 4.48e-04 157 -4.40e-04  1.33e-03   0.993  0.3224
##  lung       -2.67e-03 5.71e-04 157 -3.80e-03 -1.54e-03  -4.675  <.0001
## 
## degree = cubic:
##  BREATH O2LEVEL.trend       SE  df  lower.CL  upper.CL t.ratio p.value
##  buccal     -3.50e-07 2.31e-05 157 -4.60e-05  4.53e-05  -0.015  0.9879
##  lung        6.73e-05 2.95e-05 157  9.13e-06  1.26e-04   2.285  0.0236
## 
## Confidence level used: 0.95

newdata <- with(mullens, list(
  O2LEVEL = modelr::seq_range(O2LEVEL, n = 1000),
  BREATH = levels(BREATH)
))
mullens.grid <- mullens.glmmTMB1b %>%
  emmeans(~ O2LEVEL | BREATH, at = newdata) %>%
  as.data.frame()
mullens.grid %>%
  group_by(BREATH) %>%
  summarise(value = O2LEVEL[which.max(emmean)])

Conclusions:

• the rate of buccal breathing by buccal breathing toads declines linearly with increasing oxygen concentration
• the rate of buccal breathing by lung breathing toads as oxygen concentration increases follows a quadratic trend. Initially, as oxygen concentrations increases, the rate of buccal breathing increases up to an oxygen concentration of approximately 20 percent, after which it declines.

beta

mullens.glmmTMB2b %>%
  emtrends(specs = "BREATH", var = "O2LEVEL", max.degree = 3) %>%
  summary(infer = TRUE)

## OR
mullens.glmmTMB2b %>% emtrends(specs = "BREATH", var = "O2LEVEL", max.degree = 3, infer = c(TRUE, TRUE))

## degree = linear:
##  BREATH O2LEVEL.trend       SE  df  lower.CL  upper.CL t.ratio p.value
##  buccal     -3.46e-02 1.05e-02 157 -5.54e-02 -1.38e-02  -3.290  0.0012
##  lung       -9.07e-03 1.39e-02 157 -3.66e-02  1.85e-02  -0.651  0.5161
## 
## degree = quadratic:
##  BREATH O2LEVEL.trend       SE  df  lower.CL  upper.CL t.ratio p.value
##  buccal      1.64e-04 2.13e-04 157 -2.57e-04  5.84e-04   0.769  0.4432
##  lung       -1.65e-03 3.52e-04 157 -2.34e-03 -9.53e-04  -4.683  <.0001
## 
## degree = cubic:
##  BREATH O2LEVEL.trend       SE  df  lower.CL  upper.CL t.ratio p.value
##  buccal      2.99e-06 1.46e-05 157 -2.58e-05  3.18e-05   0.205  0.8378
##  lung        4.78e-05 2.10e-05 157  6.32e-06  8.93e-05   2.276  0.0242
## 
## Confidence level used: 0.95

## newdata <- with(mullens,  list(O2LEVEL=seq(min(O2LEVEL),  max(O2LEVEL),  len=100),
##                                BREATH=levels(BREATH)))
mullens.grid <- with(mullens, list(
  O2LEVEL = modelr::seq_range(O2LEVEL, n = 1000),
  BREATH = levels(BREATH)
))
newdata <- emmeans(mullens.glmmTMB2b, ~ O2LEVEL | BREATH, at = mullens.grid) %>% as.data.frame()
newdata %>%
  group_by(BREATH) %>%
  summarise(value = O2LEVEL[which.max(emmean)])

## r.squaredGLMM(mullens.glmmTMB2a)
performance::r2_nakagawa(mullens.glmmTMB2b)

## # R2 for Mixed Models
## 
##   Conditional R2: 0.809
##      Marginal R2: 0.334

Conclusions:

• the rate of buccal breathing by buccal breathing toads declines linearly with increasing oxygen concentration
• the rate of buccal breathing by lung breathing toads as oxygen concentration increases follows a quadratic trend. Initially, as oxygen concentrations increases, the rate of buccal breathing increases up to an oxygen concentration of approximately 20 percent, after which it declines.

Summary figures

sqrt (Gaussian) with lmer

mullens.grid <- with(
  mullens,
  list(
    BREATH = levels(BREATH),
    O2LEVEL = modelr::seq_range(O2LEVEL, n = 100)
  )
)
newdata <- mullens.lmer1a %>%
  emmeans(~ O2LEVEL | BREATH, at = mullens.grid, type = "response") %>%
  as.data.frame() %>%
  mutate(across(c(emmean, lower.CL, upper.CL), function(x) x^2))
## mutate(emmean = emmean^2,
##        lower.CL=lower.CL^2,
##        upper.CL=upper.CL^2)
head(newdata)

ggplot() +
  geom_ribbon(
    data = newdata,
    aes(
      ymin = lower.CL, ymax = upper.CL,
      x = O2LEVEL, fill = BREATH
    ), alpha = 0.3
  ) +
  geom_line(
    data = newdata,
    aes(y = emmean, x = O2LEVEL, color = BREATH)
  ) +
  theme_classic()

sqrt (Gaussian)

mullens.grid <- with(
  mullens,
  list(
    BREATH = levels(BREATH),
    O2LEVEL = modelr::seq_range(O2LEVEL, n = 100)
  )
)
newdata <- mullens.glmmTMB1b %>%
  emmeans(~ O2LEVEL | BREATH, at = mullens.grid, type = "response") %>%
  as.data.frame() %>%
  mutate(across(c(emmean, lower.CL, upper.CL), function(x) x^2))
## mutate(emmean = emmean^2,
##        lower.CL=lower.CL^2,
##        upper.CL=upper.CL^2)
head(newdata)

ggplot() +
  geom_ribbon(
    data = newdata,
    aes(
      ymin = lower.CL, ymax = upper.CL,
      x = O2LEVEL, fill = BREATH
    ), alpha = 0.3
  ) +
  geom_line(
    data = newdata,
    aes(y = emmean, x = O2LEVEL, color = BREATH)
  ) +
  theme_classic()

beta

mullens.grid <- with(
  mullens,
  list(
    BREATH = levels(BREATH),
    O2LEVEL = modelr::seq_range(O2LEVEL, n = 100)
  )
)
newdata <- mullens.glmmTMB2b %>%
  emmeans(~ O2LEVEL | BREATH, at = mullens.grid, type = "response") %>%
  as.data.frame()
head(newdata)

ggplot() +
  geom_ribbon(
    data = newdata,
    aes(
      ymin = lower.CL, ymax = upper.CL,
      x = O2LEVEL, fill = BREATH
    ), alpha = 0.3
  ) +
  geom_line(
    data = newdata,
    aes(y = response, x = O2LEVEL, color = BREATH)
  ) +
  scale_y_continuous("Buccal breathing rate", labels = function(x) 100 * x) +
  theme_classic()

obs <- mullens %>%
  mutate(
    .fixed = predict(mullens.glmmTMB2b, re.form = NA),
    .resid = residuals(mullens.glmmTMB2b),
    PartialObs = plogis(.fixed + .resid)
  )

ggplot() +
  geom_point(data = mullens, aes(y = pzBUC, x = O2LEVEL, color = BREATH), alpha = 0.2) +
  geom_point(data = obs, aes(y = PartialObs, x = O2LEVEL, color = BREATH)) +
  geom_ribbon(
    data = newdata,
    aes(
      ymin = lower.CL, ymax = upper.CL,
      x = O2LEVEL, fill = BREATH
    ), alpha = 0.3
  ) +
  geom_line(
    data = newdata,
    aes(y = response, x = O2LEVEL, color = BREATH)
  ) +
  scale_y_continuous("Buccal breathing rate", labels = function(x) 100 * x) +
  theme_classic()

References